

This example implements the coupled Cahn-\/\-Hilliard and Allen-\/\-Cahn equations for phase-\/field modeling, as described by the following weak form of the P\-D\-E. The two scalar fields are composition, $c$, and an order parameter, $\eta$). Note the application of the higher-\/order Dirichlet boundary condition $\nabla c\cdot\boldsymbol{n}=0$ using Nitsche's method.

Cahn-\/\-Hilliard\-:

\begin{eqnarray*} 0 &=& \int_\Omega \left(w_1\frac{c - c_{prev}}{\mathrm{d}t} + M\left(\nabla w_1\cdot(f_{,cc}\nabla c + f_{,c\eta}\nabla\eta) + \kappa_1\nabla^2 w_1\nabla^2 c\right)\right) dV\\ &\phantom{=}& - \int_{\partial\Omega} \left(w_1j_n + M\kappa_1\left(\nabla^2c(\nabla w_1\cdot\boldsymbol{n}) + \nabla^2w_1(\nabla c\cdot\boldsymbol{n})\right) - \tau(\nabla w_1\cdot\boldsymbol{n})(\nabla c\cdot\boldsymbol{n})\right) dS \end{eqnarray*}

Allen-\/\-Cahn\-:

\begin{eqnarray*} 0 &=& \int_\Omega \left(w_2\frac{\eta - \eta_{prev}}{\mathrm{d}t} + L\left(w_2 f_{,\eta} + \kappa_2\nabla w_2\cdot\nabla \eta\right)\right) dV \end{eqnarray*}

Free energy density\-:

\begin{eqnarray*} f(c,\eta) = (c - 0.1)^2(c - 0.9)^2\eta + (\eta - 0.2)^2(\eta-0.8)^2 \end{eqnarray*}

To implement this model, we will specify the following through defining user functions\-: \par

\begin{DoxyItemize}
\item Initial conditions \par

\item Constitutive model (via free energy density functions) \par

\item Parameter values \par

\item Weak form of the P\-D\-Es \par

\end{DoxyItemize}

First, we include the header file declaring the required user functions. These functions will be defined in this file.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Now, we first define any optional user functions. Optional user functions have a default definition that can be redefined by the user using a function pointer. This will be done in the {\ttfamily define\-Parameters} function. The available list of optional user functions includes\-: {\ttfamily boundary\-Conditions}, {\ttfamily scalar\-Initial\-Conditions}, {\ttfamily vector\-Initial\-Conditions}, {\ttfamily load\-Step}, {\ttfamily adaptive\-Time\-Step}, and {\ttfamily project\-Fields}. In this example, we redefine only the {\ttfamily scalar\-Initial\-Conditions} function, while using the default functions for the others.

{\bfseries  The {\ttfamily scalar\-Initial\-Conditions} function }

We initialized the composition and order parameters fields to be random about 0.\-5.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  Free energy density derivative functions }

This phase-\/field implementation requires several derivatives of the chemical free energy density function $f(c,\eta) = (c - 0.1)^2(c - 0.9)^2\eta + (\eta - 0.2)^2(\eta-0.8)^2$. We define the functions computing $\partial f/\partial\eta$, $\partial^2 f/\partial c \partial c$, and $\partial^2 f/\partial\eta \partial c$ here. Note that these free energy derivative functions are used only in this file. They are not members of any class, nor will we use them to set any function pointers.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  The {\ttfamily define\-Parameters} function }

The user is required to define the {\ttfamily define\-Parameters} and {\ttfamily residual} functions. The {\ttfamily define\-Parameters} defines variables and functions in the {\ttfamily App\-Ctx} object. The {\ttfamily App\-Ctx} object is defined in the app\-Ctx.\-h file. This function is used to define any values in {\ttfamily user} that will be needed in the problem. It is also used to set any function pointers for user functions that we have redefined.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Here, we define the mesh by setting the number of elements in each direction, e.\-g. a 20x20 element mesh.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We also define the dimensions of the domain, e.\-g. a unit square.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We can define a periodic (or partially periodic) domain. The default is no periodicity in all directions. Here, we override the default and define periodicity in both the x and y directions.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We define the time step and total simulation time. We also have the options to use restart files, in which case we would set the iteration index and time at which to start. We leave these values at zero to begin a new simulation. We also have the option to output results at regular intervals (e.\-g. every 5 time steps).


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We specify the number of vector and scalar solution and projection fields by adding the name of each field to their respective vector. Here, we have two scalar solution field (the composition and an order parameter). We do not use any vector solution fields or projection fields in this example.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We can specify the polynomial order of the basis splines, as well as the global continuity. Note that the global continuity must be less than the polynomial order. Here, we use quadratic basis functions with C-\/1 global continuity.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Finally, we redirect the desired user function pointers to the {\ttfamily scalar\-Initial\-Conditions} function that we defined above. This completes the {\ttfamily define\-Parameters} function.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  The {\ttfamily residual} function }

The residual function defines the residual that is to be driven to zero. This is the central function of the code. It is set up to follow the analytical weak form of the P\-D\-E. It has a number of arguments that give problem information at the current quadrature point.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\ttfamily d\-V} is a boolean, \char`\"{}true\char`\"{} if {\ttfamily residual} is being called for the volume integral and \char`\"{}false\char`\"{} if {\ttfamily residual} is being called for the surface integral.\par
{\ttfamily d\-S} is a boolean, \char`\"{}false\char`\"{} if {\ttfamily residual} is being called for the volume integral and \char`\"{}true\char`\"{} if {\ttfamily residual} is being called for the surface integral.\par
{\ttfamily x} gives the coordinates of the quadrature point.\par
{\ttfamily normal} gives the unit normal for a surface quadrature point.\par
{\ttfamily c} gives the information (values, gradients, etc.) for the scalar solution fields at the current quadrature point (see documentation for solution\-Scalars class).\par
{\ttfamily u} gives the information (values, gradients, etc.) for the vector solution fields at the current quadrature point (see documentation for solution\-Vectors class).\par
{\ttfamily w1} gives the information for the scalar test functions.\par
{\ttfamily w2} gives the information for the vector test functions.\par
{\ttfamily user} is a structure available for parameters related to the initial boundary value problem (e.\-g. elasticity tensor).\par
{\ttfamily r} stores the scalar value of the residual for the weak form of the P\-D\-E which is then used by the core assembly functions.

The following functions are available for the solution objects {\ttfamily c} and {\ttfamily u}, where the argument is the field index, i.

{\ttfamily c.\-val(i)} -\/ Value of scalar field i, scalar \par
{\ttfamily c.\-grad(i)} -\/ Gradient of scalar field i, 1st order tensor \par
{\ttfamily c.\-hess(i)} -\/ Hessian of scalar field i, 2nd order tensor \par
{\ttfamily c.\-laplacian(i)} -\/ Laplacian of scalar field i, scalar \par
{\ttfamily c.\-val\-P(i)} -\/ Value of scalar field i at previous time step, scalar \par
{\ttfamily c.\-grad\-P(i)} -\/ Gradient of scalar field i at previous time step, 1st order tensor \par
{\ttfamily c.\-hess\-P(i)} -\/ Hessian of scalar field i at previous time step, 2nd order tensor \par
{\ttfamily c.\-laplacian\-P(i)} -\/ Laplacian of scalar field i at previous time step, scalar

{\ttfamily u.\-val(i)} -\/ Value of vector field i, 1st order tensor \par
{\ttfamily u.\-grad(i)} -\/ Gradient of vector field i, 2nd order tensor \par
{\ttfamily u.\-hess(i)} -\/ Hessian of vector field i, 3rd order tensor \par
{\ttfamily u.\-val\-P(i)} -\/ Value of vector field i at previous time step, 1st order tensor \par
{\ttfamily u.\-grad\-P(i)} -\/ Gradient of vector field i at previous time step, 2nd order tensor \par
{\ttfamily u.\-hess\-P(i)} -\/ Hessian of vector field i at previous time step, 3rd order tensor

Similar functions are available for the test functions. Also, the following tensor operations are useful\-:

Tensor operations\-: \par
{\ttfamily operator+} -\/ tensor addition \par
{\ttfamily operator-\/} -\/ tensor subraction \par
{\ttfamily operator$\ast$} -\/ single contraction between tensors or scalar multiplication \par
{\ttfamily double\-\_\-contract} -\/ double contraction of two 2nd order tensors, or a 4th order tensor and a 2nd order tensor. \par
{\ttfamily trans( )} -\/ transpose 2nd order tensor \par
{\ttfamily trace( )} -\/ trace of 2nd order tensor \par
{\ttfamily det( )} -\/ determinant of 2nd order tensor \par
{\ttfamily inv( )} -\/ inverse of 2nd order tensor \par
 The example code here implements the weak form for the coupled Cahn-\/\-Hilliard and Allen-\/\-Cahn equations, as shown above.

First, we set the values for necessary parameters.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Next, we get the values for the free energy derivatives $f_{,cc}$, $f_{,c\eta}$, and $f_{,\eta}$ based on the current quadrature point.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Now, we compute the residual in a manner very similar to the analytical form. First, the weak form of the Cahn-\/\-Hilliard equation\-:

\begin{eqnarray*} 0 &=& \int_\Omega \left(w_1\frac{c - c_{prev}}{\mathrm{d}t} + M\left(\nabla w_1\cdot(f_{,cc}\nabla c + f_{,c\eta}\nabla\eta) + \kappa_1\nabla^2 w_1\nabla^2 c\right)\right) dV\\ &\phantom{=}& - \int_{\partial\Omega} \left(w_1j_n + M\kappa_1\left(\nabla^2c(\nabla w_1\cdot\boldsymbol{n}) + \nabla^2w_1(\nabla c\cdot\boldsymbol{n})\right) - \tau(\nabla w_1\cdot\boldsymbol{n})(\nabla c\cdot\boldsymbol{n})\right) dS \end{eqnarray*}


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


And finally, the Allen-\/\-Cahn equation, which completes the residual function.

\begin{eqnarray*} 0 &=& \int_\Omega \left(w_2\frac{\eta - \eta_{prev}}{\mathrm{d}t} + L\left(w_2 f_{,\eta} + \kappa_2\nabla w_2\cdot\nabla \eta\right)\right) dV \end{eqnarray*} 
\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Finally, we include a file that instatiates the template functions {\ttfamily define\-Parameters} and {\ttfamily residual}. This bit of code will generally be the same for any problem (unless you decide to use a different automatic differentation library), the user does not need to modify it.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


\section*{The complete code }


\begin{DoxyCodeInclude}
\end{DoxyCodeInclude}
 